#!/usr/bin/python -tti


######### standard library imports ###########################################

import sys
import logging
import logging.config
import os
import numpy as np

from optparse import OptionParser


######### project library imports ############################################
import settings

import readcorsika
import absorb_photons
import reflector   

import turnIntoZenit
import photonDetection

scriptpath = os.path.dirname(os.path.realpath(__file__))
logging.config.fileConfig(scriptpath+'/logging.conf')
logger = logging.getLogger('simulation_root')


parser = OptionParser()
parser.add_option("-i", "--inpath", dest="filename",
    help="MMCS file to read from", metavar="FILE")
parser.add_option("-o", "--outpath", dest="outfilename",
    help="FITS file to write into", metavar="FILE")
parser.add_option("-p", "--plot",
    action="store_true", dest="create_plots", default=False,
    help="create some informative debugging plots")
options, args = parser.parse_args()
if options.filename is None:
    logger.fatal('Please submit a filename. %s -h for help' % (sys.argv[0]))
    sys.exit(1)
if options.outfilename is None:
    logger.warning("no outpath given, no output will be generated")


##############################################################################
############# Here everything starts #########################################

thedict = {}

thedict = settings.put_settings_into_dict( thedict )

#------------ start: read the corsika file

thedict = readcorsika.readCorsikaFile( thedict, options.filename )

#------------ absorb / delete some photons

thedict = absorb_photons._absorb_photons( thedict, 
        dictname='photonData',
        absorption_key='wavelength',
        absorption_func_key = 'survival_prob_wavelength')

thedict = turnIntoZenit.turn_photons_into_zenith( thedict, 
        dictname='photonData', 
        pos_key='ground_dir', 
        dir_key='ground_pos',
        new_pos_key='ground_pos_turned', 
        new_dir_key='ground_dir_turned')

thedict = absorb_photons._absorb_photons( thedict, 
        dictname='photonData',
        absorption_key='ground_pos_turned',
        absorption_func_key = 'is_in_inner_circle_checker')

#------------ reflector simulation 

thedict = reflector.findMirrorIntersection( thedict, 
        dictname='photonData', 
        pos_key='ground_pos_turned', 
        dir_key='ground_dir_turned',
        mirror_key='mirrors',
        new_pos_key='pos_on_mirror', 
        new_mirror_id_key='mirror_id', 
        new_time_key='time_to_mirror')

thedict = reflector.reflectOnMirrors( thedict, 
        dictname='photonData', 
        mirror_alignment_point_key='mirror_alignment_point',
        pos_key='pos_on_mirror', 
        dir_key='ground_dir_turned',
        new_angle_key='angle_on_mirror', 
        new_dir_key='dir_on_mirror' )
    

thedict = reflector.findFocalPlaneIntersection( thedict, 
        dictname='photonData', 
        focal_plane_key='focal_plane',
        pos_key='pos_on_mirror', 
        dir_key='dir_on_mirror',
        new_pos_key='pos_in_focal_plane', 
        new_time_key='time_to_focal_plane')
        
thedict = reflector.findPixelChid( thedict, 
        dictname='photonData', 
        focal_plane_key='focal_plane',
        pos_key='pos_in_focal_plane', 
        new_chid_key='chid', 
        new_Coordinator_key='Coordinator', 
        new_pixel_center_key='pixel_center')
        
# At this point I guess I have to work on each event on its own...
# So let me first get a list of all event ids. I guess there can already be 
# some events missing, since some events can have zero photons in the camera.

all_event_ids = set(thedict['photonData']['event_id'].tolist())
pd = thedict['photonData']
#all_event_ids = np.array(list(all_event_ids))

if options.outfilename is not None:
    from fitswriter import RawDataFitsWriter2
    fitsfile = RawDataFitsWriter2(options.outfilename, thedict)
else:
    fitsfile = None
    
    
PulseIntegral = thedict['settings']['SipmInfo']['PulseIntegral']    # in mV s
MeanPulseWeight = thedict['settings']['SipmInfo']['MeanAmplitude']  # in "number of avalanches"
Rate = (thedict['settings']['PhotonNoiseInfo']['NsbRate']+thedict['settings']['PhotonNoiseInfo']['DcRate']) # in Hz

ac_coupling_effect = PulseIntegral * Rate * MeanPulseWeight     # in mV
logger.debug("AC coupling offset: %f", ac_coupling_effect)

pulse_function = thedict['settings']['SipmInfo']['standard_pulse_function']

NoiseLevel = thedict['settings']['AdcDict']['NoiseLevel']

AdcFullScale = thedict['settings']['AdcDict']['AdcFullScale']
AdcOffset = thedict['settings']['AdcDict']['AdcOffset']
AdcResolution = thedict['settings']['AdcDict']['AdcResolution']


SimulatedTime = thedict['settings']['SimInfo']['SimulatedTime']
SamplingFreqency = thedict['settings']['SimInfo']['SamplingFreqency']
StartTimeOffset = thedict['settings']['SimInfo']['StartTimeOffset']

NumberOfSamples = int(SimulatedTime*SamplingFreqency)
StartTimeOffsetSamples = int(StartTimeOffset*SamplingFreqency)

    
    
    
for event_id in all_event_ids:
    logger.debug("event_id %d" % (event_id,))
    
    ids = np.where(pd['event_id'] == event_id)[0]
    chids = pd['chid'][ids].astype(int)
    times = pd['time'][ids]-pd['time_to_mirror'][ids]+pd['time_to_focal_plane'][ids]
    if len(times):
        times -= times.min()
    ## The times are in ns, so at this point we convert them to sec.
    times *= 1e-9


    # Now I make a numpy array, which should hold all arrival times of all 
    # photons regardless of their origin
    # Since we only simulate a about 300ns of time with a photon rate of about
    # 100MHz lets say we need only 200 spaces for arrival times
    arrivalTimes = np.empty( (1440, 200), dtype=np.float32 ) * np.nan
    # we also want to store the weight for each photon
    weights = np.empty( (1440, 200), dtype=np.float32 ) * np.nan
    # and we store the origin as an integer code, lets say:
    #   0 : undefined
    #   1 : cherenkov
    #   2 : NSB
    #   3 : Dark Count
    #   4 : general Noise :-)
    origins = np.zeros( (1440, 200), dtype=np.uint8 )
    
    # Okay now we put the times into the arrival time containers
    for chid in range(1440):
        this_chids_times = times[chids==chid]
        N = len(this_chids_times)
        arrivalTimes[chid,0:N] = this_chids_times
        origins[chid,0:N] = 1   # <--- we handle cherenkov photons currently
    
    # Okay now we will add Nsb and DarkNoise
    NsbRate = thedict['settings']['PhotonNoiseInfo']['NsbRate']     # stored in Hz
    DcRate = thedict['settings']['PhotonNoiseInfo']['DcRate']       # stored in Hz
    Rate = (NsbRate + DcRate)       
    Period = 1./Rate                                                # in sec
    SimulatedTime = thedict['settings']['SimInfo']['SimulatedTime'] # in sec
    numberOfNoisePhotons = int(SimulatedTime * Rate)
    for chid in range(1440):
        timesOfNoisePhotons = np.random.exponential(Period, size=5*numberOfNoisePhotons)
        timesOfNoisePhotons = timesOfNoisePhotons.cumsum() - SimulatedTime
        timesOfNoisePhotons = timesOfNoisePhotons[timesOfNoisePhotons <= SimulatedTime]
        start = np.where(~np.isnan(arrivalTimes[chid]))[0]
        if len(start) > 0:
            start = start[0]
        else:
            start = 0
        N = len(timesOfNoisePhotons)
        arrivalTimes[chid,start:start+N] = timesOfNoisePhotons
        origins[chid,start:start+N] = 4
        
    # Okay and now we simulate the weight of each photon, which can be 
    # larger than 1, due to optical crosstalk
    weight_generator = thedict['settings']['SipmInfo']['AmplitudeGenerationFunction']
    weights = weight_generator(thedict, size=(1440,200))

    # Okay and now we have to think about the trigger:
    coincidenceWindow = 13*1e-9
    Threshold=20
    Trigger=False

    for patch in range(160):
        patchTimes = None
        patchWeights = None
        for pixel in range(9):
            chid = patch*9+pixel
    
            notnan_ids = np.where(~np.isnan(arrivalTimes[chid]))[0]
            arrivals = arrivalTimes[chid][notnan_ids]
            departures = arrivalTimes[chid][notnan_ids] + coincidenceWindow
            arrivalWeights = weights[chid][notnan_ids]
            departureWeights = weights[chid][notnan_ids] * -1
            alltimes = np.concatenate( (arrivals,departures) )
            allWeights = np.concatenate( (arrivalWeights,departureWeights) )
            if patchTimes is None:
                patchTimes = alltimes
                patchWeights = allWeights
            else:
                patchTimes = np.concatenate( (patchTimes,alltimes) )
                patchWeights = np.concatenate( (patchWeights,allWeights) )
                
            
        sorter_ids = np.argsort( patchTimes )
        patchWeights = patchWeights[sorter_ids]
        TriggerPathLevel = patchWeights.cumsum()
        if (TriggerPathLevel > Threshold).any():
            Trigger=True
            break
    
    logger.debug("   Trigger %s" % (str(Trigger),))
##############################################################################
    if Trigger:
        logger.debug("MakeAndSaveAdcData starting up")
    
        
        rawdata = np.zeros( (1440, NumberOfSamples) )
        x = np.linspace(StartTimeOffset,SimulatedTime-StartTimeOffset,NumberOfSamples,False)
        
        for chid in range(1440):
            notnan_ids = np.where(~np.isnan(arrivalTimes[chid]))[0]
            Ts = arrivalTimes[chid][notnan_ids]
            Ws = weights[chid][notnan_ids]
            for t,w in zip(Ts,Ws):
                pulse_function(t,rawdata[chid],w)
                
    
        ##################################### AC coupling ########################
        rawdata -= ac_coupling_effect
        
        ############### Electronics noise ########################################
        rawdata += np.random.normal( loc=0., scale=NoiseLevel, size=rawdata.shape)
    
        ############ ADC simulation, i.e. digitization ###########################        
        adcdata = rawdata.copy()
        adcdata += AdcOffset
        adcdata *= AdcResolution/AdcFullScale
        adcdata = adcdata.astype('int16')
            
        ############### Save event to fitsfile ###################################
            # Format
            #~ ('EventNum', '<u4'), 
            #~ ('TriggerNum', '<u4'), 
            #~ ('TriggerType', '<u4'), 
            #~ ('NumBoards', '<u4'), 
            #~ ('Errors', '|u1', 4), 
            #~ ('SoftTrig', '<u4'), 
            #~ ('UnixTimeUTC', '<u4', 2), 
            #~ ('BoardTime', '<u4', 40), 
            #~ ('StartCellData', '<u2', 1440), 
            #~ ('StartCellTimeMarker', '<u2', 160), 
            #~ ('Data', '<i2', (1440,1024))]) )
        if fitsfile is not None:
            fitsfile.rawdata['Data'] = adcdata
            fitsfile.rawdata['EventNum'] = event_id
            fitsfile.rawdata['NumBoards'] = 40
            fitsfile.write()
    else:
        pass

